Extended Recursion in Operator Space (EROS), a new impurity solver for the single 

impurity Anderson model 



Jean-Pierre Julien 

Theoretical Division, Los Alamos National Laboratory, MS B262, Los Alamos, NM 87545, USA and 
Institut Neel CNRS and Universite J. Fourier 25 Avenue des Martyrs, BP 166, F-38042 Grenoble Cedex 9, France 



oo 

o 

O 

(N 
> 

o 

o 

o 



I 

o 
o 



> 

(N 
O 
cn 
m 



oo 

o 



R. C. Albers 

Theoretical Division, Los Alamos National Laboratory, MS B262, Los Alamos, NM 87545, USA 

(Dated: October 20, 2008) 

We have developed a new efficient and accurate impurity solver for the single impurity Anderson 
model (SIAM), which is based on a non-perturbative recursion technique in a space of operators and 
involves expanding the self-energy as a continued fraction. The method has no special occupation 
number or temperature restrictions; the only approximation is the number of levels of the continued 
fraction retained in the expansion. We also show how this approach can be used as a new approach 
to Dynamical Mean Field Theory (DMTF) and illustrate this with the Hubbard model. The three 
lowest orders of recursion give the Hartree-Fock, Hubbard I, and Hubbard HI approximations. A 
higher level of recursion is able to reproduce the expected 3-peak structure in the spectral function 
and Fermi liquid behavior. 

PACS numbers: 71.10.-l-x, 71.20.Cf, 64.60.Cn 



In the last decade, Dynamical Mean Field Theory 
(DMFT) [H, 3 has become a widely used method to 
study strongly correlated electrons systems. It can be 
formulated by an action formalism with a self-consistent 
mapping onto a single impurity Anderson model (SIAM) . 
Because the critical step in this method is the quality 
of the impurity solvers, there has recently been an in- 
creased interest in new and more accurate SIAM solvers. 
Currently available solvers include: iterated perturba- 
tion theory (IPT) Q, which works well for small U and 
single band, the non-crossing approximation(NCA) [3], 
which is able to give the coherent peak but fails to repro- 
duce Fermi liquid behavior, and the equation of motion 
(EOM) method 0, which requires a decoupling scheme 
and misses the Kondo peak at finite U for particle-hole 
symmetric systems. IPT also presents pathologies away 
from half-filling. A numerical renormalization group 
(NRG)t0| has also been been employed but has the lim- 
itation that it captures correctly the low energy physics 
but with a logarithmic divergence instead of the expected 
Lorentzian shape of a Fermi liquid. It is also inaccurate 
for high energy physics and is limited at T=0 though fi- 
nite T is tried to be implemented Quantum Monte 
Carlo (QMC) Q, as a non-perturbative approach, is in 
principle rigorous and provides good results at high tem- 
perature for static properties but introduces uncertainty 
for dynamical properties (like spectral density) because 
of the poorly defined inverse problem, and cannot reach 
zero temperature. 

The purpose of an Extended Recursion method in Op- 
erator Space (EROS) is to provide an efficient (rapid) and 
accurate method to solve the quantum impurity problem. 
Efficiency is especially important for LDA-I-DMFT cal- 
culations, since most of the computational time is spent 



on the quantum impurity problem. An accurate solu- 
tion must include a reliable calculation for both low and 
high temperatures, low and high-energy physics, as well 
as for any filling factor for the correlated bands. We also 
intend to extend this method to quantum transport prob- 
lems for low-dimensional correlated systems (e.g., quan- 
tum dots and molecular electronics). As shown below, 
EROS is able to retrieve well-known approximations such 
as Hartree-Fock, Hubbard I, and Hubbard HI, at the low- 
est level of recursion and the Fermi liquid regime and 
metal-insulator transition at higher levels. 

The solution involves calculating a retarded Green 
fimction G{oj) =<< A; B >>ui, which is the Fourier 
transform of G{t) = ~i0{t) < {A{t),B^O)} > for the 
operators A and B. This can be done through a recur- 
sion process (RP) (see @,[ii), which can be directly seen 
by examining the coupled equations of motion for the re- 
tarded Green function for a Hamiltonian H: 



« A;B »^ 



< {A,B^ > 



<{[A,H],B^> , <{[[A,H],H],B^> 



(1) 



This can be identified with a moment expansion through 
the definition 



=<{[... [A, St} >, 



(2) 



where H appearing n times introduces a sttper-operator 
7i, the 'Liouvillian', acting in operator space: 



AH = [A, H] 



(3) 



such that /.in =< {AW\B^} >. From its moment ex- 
pansion it is possible to reconstruct the Green's function 
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as a continued fraction (CF) [l^], but a better condi- 
tioned approach is to directly obtain the CF coefficients 
from the recursion method of Haydock ll|, which can 
be directly applied in the space of operators, with the 
Liouvillian playing the role of the Hamiltonian in the 
standard RP. The more conventional wave-function RP 
consists in starting from an initial chosen state ipQ and 
then generating a set of orthonormal states ipn {n > 0) 
by the recurrence relation (60 = 0): 



(4) 



The coefficients a„ and hn+i are obtained by the scalar 
product of tpn by Htpn and with the norm || Hipn—a,ni^n — 
bnipn-i II respectively. In the ■(/)„ basis, the Hamiltonian 
has a tridiagonal form, which can be represented as a 
semi-infinite tight-binding (TB) chain. Once a sufficient 
number of coefficients pairs have been calculated, the di- 
agonal clement of the resolvent {ui — H)~^ can be ex- 
pressed as a CF. Standard techniques exist to terminate 
this expansion [isj . The application of this approach to 
operator space is straightforward, each vector now corre- 
sponding with an operator. For our purposes we will be 
mainly concerned with the local Green function, where 
A = B = cq, the destruction operator of an electron 
of given spin at site 0. The action of the Liouvillian 
on any creation or destruction operator (or any combi- 
nation of them) is then easily computable from its def- 
inition, Eq. ([3]). A natural choice for the scalar prod- 
uct, which is necessary to compute the CF coefficients, 
is {A\B) =< {A, fit} > as suggested by Eq. H]), where 
thermal averages reduce to their ground-state expecta- 
tion value at zero temperature. Details for computing 
them will be given in a forthcoming paper. In this for- 
malism, the Green's function appears as the diagonal el- 
ement of the resolvent of the Liouvillian: 

G{lo) ~« Co; Co (co(cj - H)"^|co) (5) 

This expression is the origin of the name we have pro- 
posed for this approach: the Extended Recursion in Op- 
erator Space (EROS); it reduces to the usual RP when 
only one-body operators appear in the Hamiltonian. 
We now apply it to the SIAM Hamiltonian: 

HsiAM = eo ^ "-Oo- + U notnox 

(T 
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which describes a localized orbital that has electronic 
correlations because of the Coulomb interaction U and is 
coupled to itinerant non-interacting electrons of disper- 
sion £fc; the latter is reflected by an hybridization function 
A (a;) = uj^e^. ' when spin indices arc ignored. If wc 
use a conventional RP to tridiagonalizc the last 2 terms, 
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FIG. 1; Schematic representation of the Liouvillian 



we can rewrite Eq. ([6]) as: 

H""P = Uno^noi + ^ £0^10.7 

+ {o^pCpa'^pa + f3pCp_i^Cpa + Pp+lCp^i^Cpa) (7) 

p>0o- 

These coefficients {ap,/3p+i}, are the tight-binding pa- 
rameters of a semi-infinite chain and also those of the 
CF expansion considered as a parametrization of the hy- 
bridization function A(a;), as already often noticed in 
NRG context [l3|)j whereas a RP applied to Hamilto- 
nian ([7]) provides the CF coefficients {a„,6„+i} of the 
impurity Green function G{uj). Eqs. ^ and ([7]) give the 
action of Ti. on the basis operators that is necessary to 
perform an operator RP (where the up-spin is assumed 
if a spin- index is omitted): 



coH = eoco + Pici + UcqcI^Cq^ 



+ [ipCp-i + (5p+iCp+i p > 0. 



(8) 
(9) 



Compared to the "natural" operators Cp, related to a 
given state p, the local 2-body interaction T-Lmt due to 
UriQ^nQi has generated a new operator CgCQ|Cg|, which 
is the origin of a site representation for 1/8 of a sim- 
ple tri-dimensional cubic lattice, where each site repre- 
sents an operator CpC^^^^c^^ indexed by 3 positive integers 
p,q,r (see Fig. [1}. This construction is supported by 
the following observation, which can be systematically 
extended: 

coUqiH = UcoUqi + /3iCicJ^Co| - /3iCqc[^Cq^ 

+/3iCo4xCix (10) 

The RP can be performed in this landscape as easily as 
the regular recursion on an usual lattice through knowing 
the action of 7i on operators like CpcJ^c^^. Special atten- 
tion should be paid for the case where one or two of the 
indexes p,q,r arc zero, viz., the 3 faces and the 3 edges 
of the 1/8 of the cubic lattice, since in these cases Hint 
has a non-zero effect. For the edges, this operation gives 
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b) G(co)= 0^10^!0^ig^ 

FIG. 2: Two diflerent recursion processes for G(cj) 

either or [/ for the on-site term, but for the faces, a 
new kind of operator is generated, which is the product 
of 5 "natural" operators (3 destruction and 2 creation 
operators); and the action of the LiouviUian (which are 
not representable by Fig. [1]) can be represented as the 
sites of 5-dimensional hypercube. These operators corre- 
spond to the motion of a hole with 2 electron-hole pairs. 
In our calculations we did not take them exactly into ac- 
count, but instead projected them onto existing products 
of 3 operators (which corresponds to a certain EOM de- 
coupling) or simply neglected them; in term of moments 
their influence only starts with the 7th moment. 

The self-energy, S(a;), which is related to the Green's 
function through 

G(w) = (tj - S(w) - A(w))-\ (11) 

can be expanded as a CF: 

I](c.) = A, + . (12) 

uo — Ai • • • 

To determine the coefficients of this CF, we note that the 
right hand side of Eq. (jlip , which includes the sum of A 
and E in the denominator, can be expanded in a CF by 
a RP on their 2 representative chains coupled by their 
common first site as shown in Fig.[2l[a) [l3|; at the same 
time this also has to give the the CF expansion for G in 
Fig. [UJb), which is calculated from the RP in operator 
space. This enables the coefficients of Eq. ((T2|) for the 
self-energy E to be determined. 

We now describe an application of our EROS method- 
ology to a DMFT solution of the Hubbard model. An 
excellent description of the DMFT method is given in 
[l[ and a recent review Q describes recent developments 
and extensions. The Hubbard model on a lattice (for site 
indices i and j) can be written as 

H tijc\^Cja + UY^ iM]n,i (13) 

The motivation for DMFT arises from the observation 
based on a diagrammatic analysis that in infinite di- 
mensions [l3| the self-energy becomes local, i.e., k- 
independent or simply E(a;). In lower dimensions, this is 
often only approximately true. In the local limit for E the 



Hubbard model reduces to the problem of a correlated 
impurity at site "0" embedded in an effective medium 
where all the effects of correlation are represented by the 
self-energy E(w). This can be considered as a complex 
and energy-dependent on-site energy in terms of a strong 
analogy with the TB language used in the Coherent Po- 
tential Approximation (CPA, for a detailed discussion see 
[lit), which was developed for studying alloys and has 
now been further extended for strongly correlated elec- 
tron systems [l^ : 

HdMFT = C/?l0T"-0i+ ^ Ujc\^Cja-^ ^ S(cj)71ic, (14) 

It can be shown that this problem maps onto a SIAM 
model, Eq. ([6]). To understand this within our approach, 
we start by noticing that the second term of Eq. ([T4|) , the 
"hopping" term, can be tridiagonalized by a RP to pro- 
vide the CF expansion of the "bare" hybridization func- 
tion A(tt'). The last term then functions like a constant 
onsite energy in a TB Hamiltonian that simply shifts the 
frequency by S(<x'). Hence the last two terms of Eq. (fT4|) 
can be represented by a Green's function for an effective 
TB model with the bare hybridization function A(a;) re- 
placed with an effective A(w) = A(w — E(w)). To com- 
plete the DMFT method, we then require that this self- 
consistent condition be supplemented by the requirement 
that the local impurity Green's function, computed from 
the Hamiltonian of Eq. p^. be equal to the effective 
lattice Green's function, which has E(a;) on all sites, in- 
cluding site "0". In RP language, this causes the the 
local E(w) to be added to the on-site terms ap in the 
chain representation for A(ti;). Including this self-energy 
at site "p" is equivalent to attaching to the TB semi- 
infinite chain the CF expansion of S as parameters (see 
Appendix of [l3|), and leads to a comb-shape topology for 
computing A(w). Having performed again a RP on this 
object provides a CF expansion for A(ti;). The DMFT 
self-consistency is achieved by simultaneously satisfying 
all of these CF equations as one steps through the recur- 
sion procedure. The first few recursion steps generate the 
CF coefficients Aq = U< no-a > of E, corresponding to 
the Hartree-Fock (HE) energy-independent self-energy, 
and then Bi = Un(n — 1) and Ai ~ U{1 — n) gives 
Hubbard I, which is exact in the atomic limit. If only 
operators CpcJ^coi on one edge of the cubic lattice are 
retained, Hubbard HI is recovered. Further recursions 
generate approximations that go well beyond these (see 
below). It is worth mentioning that our approach has 
some similiarities with the approach developed in [20| . 
but our direct use of a RP enables us to go beyond this 
work by avoiding their requirement to orthogonalize each 
new operator, which becomes becomes difficult after the 
first few steps. 

We have benchmarked our method on the well-studied 
half-filled case. For simplicity, all calculations were per- 
formed at zero temperature and used as a HF solution of 
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FIG. 3: RP DMFT calculation of the Hubbard model for the 
spectral density vs energy for different values of U (color on- 
line) 

Eq. (|14p as an approximate ground state for computing 
the scalar products of the operators. In the future we will 
try to use a better ground state such as a Gutzwiller ap- 
proximation or an exact diagonalization for a given num- 
ber of sites. By including approximate excited states in a 
thermal average, finite-temperature calculations are also 
possible. For the lattice model that describes the itin- 
erant electrons, we used a semi-elliptic non-interacting 
density of states, with the energy scale set by the band- 
width. Other lattices would only change the input CF 
coefficients for the hybridization function. For several 
different values of U, the spectral density is displayed 
in Fig. [3l One clearly observes the characteristic three- 
peak features in the Fermi liquid regime: the coherent 
contribution of the central peak, and the lower and up- 
per Hubbard subbands with a splitting of the order of 
U . One also observes a redistribution of the spectral 
weight of the coherent central peak as the interaction U 
increases. The Fermi liquid behavior and the intermedi- 
ate regime can be better monitored by considering the 
real and imaginary parts of the self-energy. In the Fermi 
liquid regime, we have checked that the imaginary part 
of the self-energy vanishes as it should be at the Fermi 
level Ep, and that it has the usual {uj — EpY behavior 
around Ep. 

We have presented a new impurity self-consistent 
solver for the SIAM model and a DMFT solution of 
the Hubbard model using a RP procedure in an op- 
erator space. The lowest order approximations to this 
method generate the Hartree-Fock, Hubbard I, and Hub- 
bard HI approximations. Increasing the continued frac- 
tions go well beyond these solutions, generating the co- 
herent Kondo peak, for example, and can provide increas- 
ingly accurate results. The way that self-consistency is 
achieved, step by step, makes the method efficient and 
very promising for its use in DMFT approaches and other 
more complex correlation problems. 

JPJ is grateful to J.X. Zhu (Los Alamos) for helpful 
discussions, and F. Harris (Gainesville) for sharing with 
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